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Abstract. While the cosmic baryon asymmetry has been measured at high accuracy to be 
6.1 X 10^^", a corresponding lepton asymmetry could be as large as 10"'^ if it hides in the 
neutrino sector. It has been known for some time that such an asymmetry could be generated 
from a small initial asymmetry given the existence of a sterile neutrino with a mass less than 
the mass of the active neutrino. While the magnitude of the final lepton asymmetry is 
deterministic, its sign has been conjectured to be chaotic in nature. This has been proven 
in the single momentum approximation, also known as the quantum rate equations, but 
has up to now not been established using the full momentum dependent quantum kinetic 
equations. Here we investigate this problem by solving the quantum kinetic equations for a 
system of 1 active and 1 sterile neutrino on an adaptive grid. We show that by increasing 
the resolution, oscillations in the lepton asymmetry can be eliminated so the sign of the final 
lepton asymmetry is in fact deterministic. This paper also serves as a launch paper for the 
adaptive solver LASAGNA which is available at http://users-phys.au.dk/steen. 
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1 Introduction 

The possible existence of sterile neutrinos, new light, SU{2) x U{1) gauge singlet fermions, has 
recently received a lot of attention [1]. The reason is that an additional neutrino mass state 
seems to be needed in order to explain all available neutrino oscillation data. In particular, 
short baseline oscillation data seems to point to the existence of one or more mass states 
beyond the three predicted by the standard model, and the same is true for reactor neutrino 
flux measurements. However, any additional neutrino states must be flavour singlet states 
because of the well established LEP bound on the number of neutrino flavour states. The 
oscillation data points to a relatively small mixing angle between the active sector and the 
new sterile sector which means that the additional mass eigenstates are close to being true 
flavour singlet states. However, even a very small mixing can cause the new mass states 
to be equilibrated in the early universe by a combination of scattering and mixing [2, 3] - 
a phenomenon which will have a profound impact on cosmological observables such as the 
cosmic microwave background. For example the excess energy density carried by the new 
mass states might be observable in the CMB because of an increase in the early integrated 
Sachs- Wolfe effect [4]. 

In the present paper we wish to investigate a long-standing open question related to 
active-sterile neutrino oscillations in the early universe. It was realised quite early that a 
small pre-existing lepton asymmetry can be strongly enhanced by active-sterile oscillations 
to the extent that it might become important to sterile neutrino equilibration in the early 
universe [5-7]. Furthermore, in certain regions of parameter space the lepton asymmetry 
could become oscillatory in nature for the inverted hierachy case, and lead to a random 
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sign of the final lepton asymmetry - in other words the sign of the final lepton asymmetry 
exhibited all the usual signs of being a chaotic phenomenon [8] . 

The chaoticity of the lepton asymmetry generation has been firmly established for the 
so-called quantum rate equations where all neutrinos are assumed to be at the same mo- 
mentum such that the entire neutrino distribution passes through the MSW resonance at 
the same time [9]. However, it was not clear whether this phenomenon is truly physical or 
a consequence of the approximation scheme used. It is certainly conceivable that the lepton 
asymmetry oscillations become damped if the full momentum dependent quantum kinetic 
equations are used. Some steps to investigate this phenomenon were taken in Refs. [10, 11] 
and in the seminal paper [12] a new numerical scheme was developed for resolving very nar- 
row resonances in momentum space which previously made numerical convergence almost 
impossible to achieve. Even with such numerical tricks, the system of equations is extremely 
challenging to solve numerically. In this paper we combine the adaptive grid of [12] with a 
pair of solvers for stiff differential equations and a high performance multi-threaded sparse 
LU-decomposition . 

We convincingly demonstrate that with very high numerical resolution the lepton asym- 
metry oscillations are damped away and the sign of the lepton asymmetry becomes determin- 
istic. In other words, the chaos previously found in the system is of numerical, not physical 
origin. 

The paper is structured as follows: Section 2 contains a description of the quantum 
kinetic equations and the quantum rate equations. Section 3 discusses Lyapunov exponents 
as a tool for quantifying information loss in the system. Section 4 contains our main results 
and we provide our conclusions in Section 5. The Appendix contains a brief description of 
our numerical code which is now publicly available. 

2 Equations of motion 

The full system of oscillating and scattering neutrinos in the early universe at temperatures 
close to neutrino decoupling is in principle described by the full A^-body Hamiltonian which 
can be followed in time using the Liouville-Von Neumann equation for the density matrix. 
However, in the system of interest here, correlations introduced by collisions can safely be 
ignored and the system can be followed using the reduced 1-body density matrix(see [13] for 
an excellent treatment of this point). The system of equations for the reduced 1-particle den- 
sity matrix is a generalisation of the 1-particle Boltzmann equation, known as the quantum 
kinetic equations. 

2.1 Quantum kinetic equations 

We use the density matrix formalism to describe the oscillations between active and sterile 
neutrinos, and we parametrise the density matrix using Bloch vector components: 



where rr is the Pauli matrix, and /o = (1 + exp(p/T))~^ is the Fermi-Dirac distribution 
function with zero chemical potential. 

Since the lepton number is calculated from the asymmetry between neutrinos and anti- 
neutrinos, it is crucial to know this quantity with good numerical precision. We achieve this 




(2.1) 
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by changing to the symmetric and asymmetric variables 

P^ = Pi±Pi , i = 0,x,y,z. (2.2) 
Furthermore, we change to the more physical variables 

P±=Po^ + P± = 2^, (2.3) 
Jo 

P±=Po^-P± = 2^. (2.4) 
Jo 

With this choice of variables, the quantum kinetic equations finally have the form [12, 14]: 

P± =y,P± + r [2/±//o - P^] , (2.5) 
Pt=-V,P^, (2.6) 
P^" =-{Vo + Vi) P± - VlP;[ - DPt, (2.7) 

P^ = iVo + Vi) Pt + VlP^ - [Pt - Pf) - DPy^. (2.8) 

Defining the co-moving momentum x = p/T and the degeneracy parameter ^ = /u/T, are 
given by 

•^et = 1 + e^-« ^ l + e^+«' ^^"^^ 

The potentials are given by 

Vx=^sin2e„ (2.10) 

^o = -^cos20„ (2.11) 

yajV^^^^^S^ (2.13) 

where the x subscript on V denotes the x-direction in the Bloch space, not to be confused 
with the comoving momentum x = p/T. The last two potentials depend on which active 
neutrino flavour we are considering. For an electron neutrino ge = 1 + 4sec^ /{n^^ + rip^) 
while for a muon or tau neutrino g^^r = 1- All number densities n^, are normalised to 1 in 
thermal equilibrium. The effective asymmetries are given as 

i(e) = Q + 2sin2 Ow^ + Q - 2sin2 9^^ Lp - + 2L^^ + L^^ + L^^, (2.14) 

L(fj.) =L^e-j - Le - Li,^ + L^^, (2.15) 
L(t) =L(^e) - Le - L^^ + L^^, (2.16) 

where Lj = (nj — nj)Nf/N~f and Nj and are the fermion and photon number densities 
in equilibrium respectively. 
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The collisional damping can be approximated as 



2 ' 



(2.17) 



where F is the scattering rate of neutrinos with other neutrinos and electrons in the electron 
neutrino case. When evaluating the full matrix element it turns out that it can be well 
approximated by 

r = CaGlxT^. (2.18) 

The constant Ca depends on the flavour and has the values Ce ~ 1-27 and C^,t ~ 0.92. 

With all the quantum kinetic equations in place, we could in principle solve the equations 
now. However, the re-population term of equation (2.5) breaks lepton number conservation 
because of the approximate form of the scattering kernel, equation (2.18). To ensure this 
conservation, we evolve an additional differential equation for L where we have put the 
offending term to zero in the integrand by hand: 

1 



L 



(a) 



8C(3) 



dxx foV^Py 



(2.19) 



Finally, we will use the following explicit form for the degeneracy parameter, ^ = fi/T [14]: 



-27r 



sinh I — arcsinh 



18^/3C(3) 



L 



(a) 



(2.20) 
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Parametrisation of momentum space 

The quantum kinetic equations, that we want to solve, have now been presented. However, 
as it was pointed out by Kainulainen and Sorri [12], any attempt of solving the equations 
using a uniform grid in momentum space is doomed, due to the very narrow resonances that 
are present for small mixing angles. In order to resolve these features, we will use a slightly 
modified version of the parametrisation they introduced. 

As a first step we map the momentum grid, x, to the variable u such that the extremal 
point of X gets mapped to u ^ 1/2: 

^ ' ^cxt ~l~ 2^max 



u{x) = K- 



K 



(2.21) 



X -\- Xcxt 3Jmax ^lain 

This improves the resolution where the distribution has most weight, but it does not help 
much in resolving the resonances. 

The resonances can be found from the condition Vz = Vq + Vi + Vl = Q. For the inverted 
hierarchy this results in the resonances [14] 



Vi 



where A = 



For the normal hierarchy the resonances are 

Vi 

Vo 
Vi 



A - ^ A^ - 
A + - l) 



(2.22) 
(2.23) 

(2.24) 
(2.25) 
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To improve the resolution at these points, we have used the same polynomials as Kainu- 
lainen and Sorri [12]. We have however used a scheme allowing for A^res points with improved 
resolution instead of only the two they used. This was done using the parametrisation 



u{v) 



av + ai + bi {v — f ^ for v < Vc^, 

av + ai + bi{v - Vn)^ for Vci_-^ < v < v^, i = 2, . . . ,n - 1 (2.26) 

^av + an + bniv - VrJ'^ for Vc„_i < 



where Vc, = Hvr, +Vr,+^). 

The differentiability of u{v) requires that 

du(v) , du(v) 
lim — ^ = lim_ — (2.27) 

and thus b = bi = bj\/i,j. When we furthermore require continuity, and that u(0) = 0, and 
u(l) = 1, we get the equations 

ai-i + b {vc^ - = ai + b {vc, - VrJ^ , (2.28) 

= ai-6t;^^, (2.29) 

1 = an + b{l-vf.J +a. (2.30) 

This system of equations cannot be solved analytically, but using Newton's method we can 
find b and Vr^ . It is also possible to set up differential equations describing , but given that 
Newton's method converges to machine precision in just a handful of iterations, we prefer 
this method due to its increased precision. Note that speed is never an issue here, since the 
work required in this step is completely negligible compared to the work required to calculate 
the remaining part of the right hand side. 

With this parametrisation in place, the concentration of points close to the resonances 
can be increased when a is decreased. However the parametrisation introduce a complicated 
time dependence, and we have to modify QKE in order to take this into account [12]: 



f dp{T,x{T,v)) \ ^(dp\ (du\ (dv\ (dp\ 
V dT \dTj,^^\dTj^\duJ^\dvJ 



+ 1 — 1 I — I {^] . (2.31) 

T 



The first term is the ordinary QKE, while the last describes the changes induced by the 
parametrisation. The three factors in the last term can be found using different methods [12] . 
{dp/dv)T can be found using a simple stencil method. {dv/du)T = (du/dv)^^ can be found 
by differentiation of equation (2.26). Finally, (du/dT)^ can also be found by differentiating 
equation (2.26), but the result depends on dvj-JdT. To determine these, equation (2.28) and 
equation (2.30) must be differentiated, and then the resulting system of linear equations can 
be solved. Note that equation (2.29) can be used to find an expression for b and db/dT, but 
it could also be included as another equation if db/dT was included as an unknown. 

2.3 Quantum rate equations 

The QRE can be derived from QKE by applying two approximations. The first and most im- 
portant approximation is to replace the momentum by its mean value {p) = 77r^/(180(^(3))T ~ 
3.15r, assuming a thermal distribution. This greatly reduces the number of equations in the 
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Figure 1. The difference between v and x. Note that x runs from only to 1, while v covers x-values 
from to 100. 



system and gives a much simpler system to solve. Next, we neglect re-population of active 
neutrinos from the plasma. This is a good approximation provided that [15] 



\dm^\sm^{29) < 10 
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(2.32) 



This is the case for the parameters we are interested in. As a consequence Pq = 0, and 
the change from Pq and Pz to Pa and Pg does not give any advantages. Finally, L can 
be calculated directly from P^ — Pz, and since we are neglecting re-population, it is not 
necessary to evolve an independent variable to ensure lepton number conservation. In the 
next paragraphs we will sum up the result of the considerations above. 
The density matrices of the system can be written as 



(2.33) 



2' " ' 2 

Since the interesting quantity is the lepton asymmetry, we will still use the variables P^ 
Pi lb Pi to improve numerical accuracy. The time derivatives of these are: 



Pt = -(vb + yi)p± - VlP^ - DPt, 

Py = (^0 + ^i)^'.^ + VlPI - V^Pt - DP^, 



V P 



(2.34) 
(2.35) 
(2.36) 



The potentials Vq, Vi, and Vl are still given by equations (2.10)-(2.13) and the damping D 
by equation (2.17), but now x = 3.15. 
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Figure 2. Diverging orbits as described in equation (3.1). 



The lepton asymmetry can be calculated from the z-component of P: [16] 

L = 2 • - npj + 2 • Linitiai = + 2Linitiai- (2.37) 

While the above considerations give the correct equations, it might not be obvious why the 
approximations are justified. For at more rigourous derivation see e.g. Enqvist, Kainulainen, 
and Thomson [16] or McKellar and Thomson [17]. 

3 Quantifying chaos 

The stochastic behaviour of the lepton asymmetry sign has been demonstrated qualitatively 
in QRE [9, 18] and for a few points in parameter space in QKE [12]. To quantify this we 
need to use some tools from the mathematical theory of chaos as it was done for QRE in [19]. 

Defining chaos mathematically is hard; the standard definition relies on two properties. 
The first property is exponential divergence of initially nearby orbits, and the second property 
is global mixing of the system [20]. This broad definition leaves room for some non-chaotic 
systems, and still it does not cover all chaotic systems that could be said to be chaotic. 
Fortunately, our interest will be confined to the divergence of close orbits. 

The divergence of close orbits, e.g. / and g in Figure 2, is described mathematically by 
a Lyapunov exponent, A. 

+ dt) - git + dt)\ = 2^'^*|/(i) - 9it)\. (3.1) 

We use the base of 2 since we would like to measure the amount of lost information in bits. 

The amount of information in a system depends on how well the initial conditions are 
known such that more decimals correspond to more information. In this framework the 
Lyapunov exponent can be interpreted as the lost information per second. The amount of 
lost information from time ti to time t2 is thus given by 

1=1 Xdt. (3.2) 

This information loss can be compared to the amount of information we have in the initial 
conditions. Thermal fiuctuations in the early universe give a random background with AL ~ 
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10~^^ while the initial lepton asymmetry is assumed to be 10"^*^. This difference corresponds 
to about 25 bits of information [10, 19]. 

The next interesting question is how to calculate the Lyapunov exponent and the infor- 
mation loss. First of all we must generalise to more than one dimension. Given some orbit / 
we will define Vj(t) to be an orthogonal set of small perturbations. Inspired by equation (3.1) 
we define the ith Lyapunov exponent by [21, 22]: 

A.a,.*.) = ijlofe(^). (3.3) 

Actually this is the finite-time Lyapunov exponent, since we do not take the limit of At — )• oo. 
This is preferable as we do not expect the limit to give an accurate description due to the 
transient nature of the chaotic behaviour in QKE and QRE. Given this definition of the ith 
Lyapunov exponent it is straightforward to calculate the information loss: 

The naive way to calculate Vj(t) would be to take an extra point in the parameter space 
for each i that has a small initial separation to the orbit and solve the differential equations 
for all of them. This procedure fails since the system mixes and the orbits diverge. Listead we 
want to probe the local properties throughout the evolution. Given the differential equation 
y' = f(y), this can be done by considering the linearised system: 

(y + dyY « f (y) + ^ ' dy = f (y) + J • dy. (3.5) 

This leads to a differential equation governing Vi{t) that can be solved along with the other 
equations: 

v^ = J-v,. (3.6) 

The method can be improved further to give better performance [21-23], but it requires the 
calculation of all Vj. This is not feasible for QKE since it would require solving 0(64?;^^^) 
differential equations, where v^es is the number of momentum bins. Instead we found only 
the largest Lyapunov exponent by solving the QKE's along with equation (3.6) for a single 
vector vi. 

Finally there are some numerical issues to consider. The calculation of J in every 
time step requires 0{vres) evaluations of the QKE, making the calculation unfeasible again. 
Instead we make a compromise in accuracy by using the Jacobian matrix computed by the 
solver although it is not calculated as often^. Another performance issue arises regarding the 
initial value of Vj. The direction of Vj must be chosen randomly to ensure the presence of 
a leading eigenvector component. However this initial random direction gives rise to a very 
sensitive system and forces the solver to take very small steps. To circumvent this problem 
we suppress the numerical importance of the initial Vj by choosing |vj| <^ 1. Both of these 
simplifications infiuence the detailed development of I{t), but it does not change the overall 
patterns. 



^The criterion for recalculating the Jacobian is based on the convergence speed of the hnearised system. 
Thus, it is reasonable to believe that the Jacobian will be recomputed if it is too different from the correct 
Jacobian. 
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= -10-2 eV2, sm^{29) = 10"^ 
-4 I ^ ^ ^ 




equations. Sign changes are marked with dots. 
4 Results 

4.1 Quantum rate equations 

We have solved the quantum rate equations numericahy for the muon neutrino case, and 
our results are similar to what others have found [9, 18, 19]. We have furthermore followed 
the evolution of the Lyapunov vector and calculated the total information loss. An example 
that demonstrates the type of oscillations we see is found in Figure 3. Since we only follow 
one momentum state, we can predict the expected position of the oscillations by solving 
Vq + Vi + Vl = 0, and we find that this prediction holds. In Figure 3 we see a significant 
number of oscillations, but if sin^(2^) is decreased the number of oscillations decrease as well 
while an increased value of sin^(2^) results in more oscillations. The number of oscillations 
and the resulting loss of information also depends on 6m?, but in a less straight forward way. 

The amount of information that is lost for a given set of parameters is seen in Figure 4. 
On top of this, the final sign of the lepton asymmetry shown as a bright or dark shading, and 
the agreement between the two patterns is quite good. The area where the lepton asymmetry 
sign seems to be very sensitive to changes in 6m? and sin^(20) are the same areas where the 
information loss exceeds the limit for stochasticity, / = 25. The slightly high values of 
information loss that show up in the lower left corner is likely due to a loss of information 
that is not related to the lepton asymmetry since the associated Lyapunov vector, which 
indicates the variables that are responsible for the information loss, is mainly pointing in the 
Px and Py direction. 

The levels of information loss in Figure 4 show the same tendency as the Lyapunov 
analysis that was done previously [19]. The numbers are not exactly the same, but the 
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Iog(sin2(20)) 

Figure 4. The Information is shown in colour. A bright shading of the colour indicates a positive 
final lepton asymmetry and a dark shading of the colour indicates a negative final lepton asymmetry. 
All / > 25 have the same colour as this is the limit of stochasticity. 

exact results depend on the numerical scheme that is used as well as the initial temperature 
and other parameter choices. This means that the loss of information cannot be used as an 
exact gauge to distinguish stochastic and non-stochastic areas, but it can help to deepen our 
understanding of the quantum rate equations. 

4.2 Quantum kinetic equations 

The quantum kinetic equations have been solved numerically using a numerical differentiation 
formulae of order l-5(ndfl5) devised by Shampine and Reichelt [24], and again we have 
focused on the muon neutrino case. Since the system is expected to be stiff, we have reduced 
the maximal order from 5 to 2, and thereby the solver becomes L-stable. The ndf 15 method 
is implicit and uses a modified Newton's method to solve the implicit equation. This imply 
the solution of matrix equations, and due to the need of many momentum bins, they take up 
most of the computation time. To improve the performance and take advantage of multicore 
CPU's we used the SuperLU_MT package [25, 26] to solve the matrix equations. SuperLUJIT 
takes advantage of a very clever work distribution mechanism, and it depends on efficient 
BLAS-implementations for speed. It is possible to use a multithreaded BLAS implementation 
in conjunction with SuperLU, but we found it most efficient to use a sequential BLAS and let 
SuperLU_MT do all the threading. We tested three different BLAS-libraries and found Intel's 
MKL BLAS to be the fastest for our purposes. 

Our main results are shown in Figure 5-7. All three cases show that the oscillations 
present for a small number of momentum bins disappear when a sufficiently large number 
of bins are used. The first figure also confirms that the solution has converged as the three 
solutions with most bins are indistinguishable. These results are also unchanged when we 
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Figure 5. Tiie lepton asymmetry for Sm^ — — lO^^eV and sin (261) = 10^'^. Stars mark sign 
changes. Note that Vrcs = 1600, 2400, and 3200 are indistinguishable. The parameters used here are 
also used in Figure 3. 

modify the parametrisation parameters within reasonable limits, and the patterns are similar 
for all the points in parameter space that we have examined. This is contrary to the findings 
using quantum rate equations(QRE) [9, 18, 19] and to the results of two earlier studies using 
the full quantum kinetic equations(QKE) [11, 12]. The chosen oscillation parameters in 
Figure 5-7 are all in the region where QRE show many oscillations, and the parameters of 
Figure 6 are in a part of the parameter space where the final sign of the lepton asymmetry 
appears to be stochastic [9, 18]. 

The lack of oscillations can be understood from the quantum kinetic equations. Consider 
the equations for P+ and P~. 

= {Vo + V^)P+ + VlP- - \v,iP^ - P+) - DP+ (4.1) 
= -{Vo + yi)P^ - ViPy - DP+ (4.2) 
Py = (^0 + Vi)P^ + VlP^ - \v.{Pa - Ps) - DPy (4.3) 

Let us take a point in momentum space above the resonant value as it is done in Figure 8. 
Here the time derivatives are very close to zero before the resonance has passed. As the 
temperature decreases, the resonance approaches the point, and (Vq + Vi)P^ becomes less 
dominant in the equation for P+. Since — 1/2T4P+ is unchanged, this leads to a rise in P+. 
The growing value of P^ effects P^ , and P^ becomes positive as well. Finally, this affects 
P~ which becomes positive due to the term VlP^- When this initial mechanism has excited 
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Figure 6. The lepton asymmetry for 5m^ = —10 ^eV^ and sin^(26') = 10 ^. Stars mark sign 
changes. 
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Figure 7. The lepton asymmetry for 5m? = — leV^ and siT?{2d) = 10 ®. Stars mar_k sign changes. 
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Figure 8. P+ , P+ , and P" as a function ofT for x = 0.3. 5m^ = -IQ-^cV^. and sm'^{2e) = 10"'^ 
as in Figure 3 and Figure 5. 

the different parts of the density matrix, the complexity of the full equations dictates the 
detailed evolution, but eventually the damping terms will dominate, and Py' becomes zero. 

To sum up the mechanism: The amplitude of P~ depends directly on the value of L, 
and as L decreases the amplitude, P~ does the same. Since dL/dt is the integral of P~ , this 
mechanism means that as L approaches zero so does dL/dt, and there will never be a sign 
change in L. 

Despite the feedback mechanism just described, we see oscillations when the number of 
momentum bins is too small. This happens because the mechanism is a gross simplification, 
and diff'erent effects that are not considered can easily make the approximations break down. 
For the case in Figure 5 we can try to understand why the solutions show oscillation for 800 
momentum bins but no oscillations for 1600 and 3200 bins. 

In Figure 9 the evolution of dL/dt is shown as a function of temperature. As it could 
be expected, there are oscillations for 800 bins. What might be more surprising is that the 
solutions with 1600 and 3200 bins show oscillations as well, however on a much smaller scale. 
Since dL/dt is the integral of Py , we consider this variable in Figure 10. At a temperature 
of 20MeV there is no difference between the three solutions. In the first panel some small 
oscillations start to emerge, and they become dominant for the solutions of the 800 and 
1600 bin cases at a temperature of 12.5MeV. The key difference between 800 and 1600 bins 
is that the oscillations are symmetric around the correct solution in the 1600 bin case and 
asymmetric in the 800 bin case. This explains the difference in Figure 9. Notice also that the 
scale of the y-axis decreases dramatically as the temperature evolves. This is to be expected 
due to the feedback mechanism described earlier. 

When we want to pinpoint the source of these rapid oscillation further, we need to 
look at dPy /dT. This is a quite complicated quantity where the different terms cancel each 
other to a very high degree, and the oscillations could in principle originate in a complicated 
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Figure 9. The oscillations in L are also present in dL/dt as it could be expected. The parameters 
are 5m? — — 10~^eV^ and sin^(20) = 10~^ as in Figure 5, and the system has been solved for the 
same values of T . Stars mark sign changes. 
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Figure 10. 6m'^ = -lO^^eV^ 

and sin (29) = 10 ^ as in Figure 5 and 9. Py is dominated by 
prominent oscillations as the temperature falls, and this leads to the oscillations seen in Figure 9. 



interplay between the different terms. Fortunately, the T^-PjT'term seems to contribute a lot 
more to the oscillations than the other terms as it is seen in Figure 10. This means that the 
oscillations come from P~ , and as Figure 12 shows the oscillations grow gradually from some 
ill resolved features. P~ has the advantage that its derivative is very simple. It only depends 
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Figure 11. Compare this Figure to Figure 10. The oscillations in P~ come from the Vx-term in 
equation (4.3). The oscillations in this term are present even for high temperatures, but they do not 
dominate until the remaining terms shrink to 0(10"'^). 
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Figure 12. Sm"^ = -lO^^eV^ and sin^(26') = 10"'' as in Figure 5, 9, 10 and 11. There are no 
oscillations atT — 25MeV, but they emerge as the negative feature is moving away from the resonance 
as the first panel shows. These oscillations are purely numerical in origin, and they arise because the 
region below the resonance is undersampled. In the second panel the oscillations extend across the 
resonance and this gives rise to the oscillations seen in Figure 11. 



on P which has no oscillations at that high temperatures, and thus the only possibility is 
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5m'^ = -leV2, sm^{2e) = lO^^^ , wres = 800 




r[GeV] 

Figure 13. Srn^ = -lO^^gV^ 

and sin {20) = 10 ^ as in Figure 5 and Figures 9-12. Stars mark sign 
changes. The osciUations disappear as alpha rises and the negative feature from Figure 12 becomes 
better sampled. Note however that some other oscillations start to emerge for high temperatures, as 
alpha rises and the sampling of the resonance becomes worse. 

that the osciUations come from the term added to account for the parametrisation. This is 
also to be expected since an insufficient resolution in momentum space results in an inaccurate 
estimate of dP~ /dv. 

The most straightforward way to eliminate the error in dP^ /dv is to increase the number 
of momentum bins in the calculation, but another option is to change the parametrisation. 
When increasing the a parameter, less points are placed close to the resonance and more 
away from it. This is shown in Figure 13 where the convergence is achieved by increasing a 
rather than increasing v-^es- The drawback of this approach is that the resonance becomes ill 
resolved at some point, and this does also result in oscillations. A hint of this effect can be 
seen in Figure 13 for high temperatures, but it is too small to influence the final result. A 
third method to improve the sampling of the negative feature in Figure 12 is to change the 
parametrisation. This can easily be done by including a third v^^ in equation (2.26). Since 

the feature starts at the original position of the resonance, and Pr = when the resonance 
has moved away, it remains at this position at least until the lepton asymmetry starts to 
grow, and this is the period we are interested in. With this modification to the code we can 
show the convergence to no oscillations even for 5m? = — leV^ and sin^ (2(9) = 10-^ which 
is in the middle of the oscillating region found by Di Bari and Foot [11]. In general, many 
different effects may result in oscillations, and we shall not try to describe all the cases here. 
However, all the oscillations we have seen disappear when the number of momentum bins is 
high enough and all features are well resolved. 

We have tried to include a smoothing term for the sterile neutrinos to confirm the 
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oscillations seen by Kainulainen and Sorri [12]. In order to conserve the number density we 
added 



but the solution converged to our original result when was gradually reduced. It is worth 
noting that the computation time was reduced with a factor of ~ 15 when the distributions 
of the sterile neutrinos were smoothed out. Unfortunately the amount of smoothing that 
is necessary to remove the numerical noise in P^^ is so large that it significantly alters the 
behaviour of L, even when the added term conserves number density. The reason for the 
high sensitivity of L is the large cancellations that occur in several different terms, one of 
them being for dL/dt which is the integral of an oscillating function as shown in Figure 10. 

4.3 Lyapunov analysis of quantum kinetic equations 

In the previous section we have shown that the oscillating behaviour found in QRE disappear 
in QKE when an adequate number of momentum bins is used. However, it could still be 
interesting to see how this disappearance manifests itself in the Lyapunov analysis of the 
system. 

When calculating the information loss, the initial direction of the Lyapunov vector is 
chosen at random. As the computation is started, this initial direction changes, and the 
length of it grows rapidly leading to a information loss of approximately 20 bits. Since this 
information loss can not be related to a physical phenomenon, but has a purely numerical 
origin, we reset the scale of information loss accordingly. The parameters we have chosen 
to investigate are once again those of Figure 5. Since the Lyapunov analysis introduces n 
additional equations, we have chosen to vary a instead of the number of bins. The results 
are shown in Figure 14. 

The oscillating case is seen to the left, and the information loss is substantial here as 
could be expected. When a = 0.1, the information loss is associated with a Lyapunov vector 
pointing in the direction of and P^. This loss of information is obviously not related to 
an unpredictable sign of the lepton number, but it is rather related to the decoherence as the 
temperature decreases. It is interesting to note that the calculation of Lyapunov exponents 
affects the convergence of the solution. When the Lyapunov exponents were not calculated a 
a = 0.5 was needed, but when the Lyapunov exponents where included, it only took a value 
of 0.1 to converge to the non-oscillating solution. This is probably due to some slight changes 
in numerical errors while solving the equations, and it further illustrates that the oscillations 
are non-physical. 

5 Conclusions 

We have studied the evolution of the lepton asymmetry in models with active sterile neutrino 
oscillations. First, we confirmed that using the quantum rate equations and assuming the 
hierarchy to be inverted there are certain regions of the mixing parameter space where the 
lepton asymmetry can undergo sign changes and where the sign of the final asymmetry is 
chaotic in nature. 

Once the full momentum dependent quantum kinetic equations are used these oscilla- 
tions vanish as soon as the resolution is sufficiently high and the system no longer exhibits 
any signs of chaotic behaviour. We quantified the behaviour of the system using a Lyapunov 
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a = 0.01 a = 0.1 
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T T 

Figure 14. 5ni^ = — lO^^eV^ and sin^(26') = 10^^ as in Figure 5 and 9-12. The Lyapunov analysis 
has been performed for a ~ 0.01 and a = 0.1. The first plot for each a show the relative components 
of the Lyapunov vector, the second plot shows the amount of lost information, and the last plot shows 
the logarithm lepton asymmetry. 

analysis and found that there is indeed a very significant amount of information loss in the 
regions where L becomes oscillatory, to the point where the final sign does become chaotic 
in nature. Using this analysis it is also clear that this information loss is not a physical phe- 
nomenon, but rather an artifact associated with lack of numerical resolution in momentum 
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space. This settles a long standing open question in early universe neutrino phenomenology. 

A final, interesting point is that the Lyapunov analysis can also be used to study kine- 
matical decoherence of the system at temperatures below resonance. At low temperatures 
the most significant information loss is no longer associated with L, but rather with Px,Py, 
i.e. in the transverse direction associated with normal vacuum oscillations. 
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A LASAGNA 

Code available at http://users-phys.au.dk/steen. This appendix gives a short overview 
of the numerical code that we have developed. 

A.l Overview 

LASAGNA^ is a code written in C for solving stiff systems of ordinary differential equations 
(ODEs). We developed the code for solving the Quantum Kinetic Equations for active-sterile 
oscillations in the early Universe, but we stress that due to its modular structure makes, it 
is easy to implement any other set of differential equations. Input parameters are read from 
a text file'^, and output is handled by a user specified output function. This could be a 
simple function that just writes to the screen, but LASAGNA also contains custom routines for 
creating, displaying and modifying its own binary file format. The LASAGNA binary files are 
compatible with MATLAB for fast and easy visualisation of output data. 

The user can choose from 3 different ODE-solvers, ndf 15, RADAU5 and an implementa- 
tion of the embedded Runge-Kutta formula due to Dormand and Prince of order 4 and 5. 
All three solvers are plug-compatible, so it is easy to try different solvers and also possible 
to implement another solver, ndf 15 and RADAU5 are both implicit solvers, and they both use 
a modified Newton's method for solving the (possibly) non-linear system of equations. The 
implicit methods also require access to a linear algebra solver. As the most optimal linear 
algebra solver for a given problem depends on the size of the system and the sparsity of the 
Jacobian, we have implemented one dense and two sparse solvers. Like the ODE-solvers, all 
the linear algebra solvers are plug-compatible. 

A. 2 ODE-solvers 

We have included 3 ODE-solvers in LASAGNA, ndf 15, RADAU5 and an explicit Runge-Kutta 
solver which is not optimised and is primarily included for reference. The first two are both 
suitable for stiff systems, ndf 15 is a variable order, adaptive step size linear multistep solver, 
based on the Numerical Differentiation Formulae of Shampine [24] which are of order 1-5. 
Since the formulae are only L-stable for order 1 and 2, we often found it necessary to reduce 
the maximal order to 2'^. 

RADAU5 employs an implicit Runge-Kutta method of order 5 which is L-stable. It is a 
3-stage method, so the resulting system of algebraic equations are naively 3A'^ x 3N. However, 

■^Possibly an acronym for Lepton Asymmetric Sterile- Active momentum-Grid Neutrino Analyser. 
^Many thanks to Julien Lesgourgues for letting us use the parser he wrote for his Boltzmann-code CLASS. 
*This is done by setting int maxk=2 in evolver_ndf 15 . c 
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denseJJR 


sparse 


SuperLU 


Type of method 
Range of sizes 
Threaded 


Dense 
1 - 100 
No 


Sparse 
~ O(10i - 102) 
No 


Sparse 
~ O(102 - 10^) 
Yes (good scahng up to 16 cores) 



Table 1. The imphcit ODE-solvers can utihse different linear algebra solvers through a common 
interface. 

by doing clever transformations on the Runge-Kutta matrix^, Hairer and Wanner showed [27] 
that the equations separate to a x system and a 2N x 2N system. Moreover, the 2N x 2N 
system are in fact equivalent to a x system of complex numbers. The amount of linear 
algebra work in each step of RADAU5 is then just 5 times the work in each step of ndf 15. 
Which one is better will depend on size, stiffness and tolerance parameters. 

A. 3 Linear algebra solvers 

A good linear algebra solver is very important for implicit solvers, but the best solver is 
problem dependent as we have shown in Table 1. It is important to have an interface which 
is general enough to accommodate a range of solvers, but without sacrificing performance. A 
wrapper to a linear algebra solver in LASAGNA consists of the following 4 subroutines which 
are passed to the evolvers: 

int linalg_initialise (MultiMatrix *A, 

void **linalg_workspace , 
ErrorMsg error_message) ; 
int linalg_f inalise (void *linalg_workspace , 

ErrorMsg error_message) ; 
int linalg_factorise (MultiMatrix *A, 

EvolverOptions *options, 
void *linalg_workspace , 
ErrorMsg error_message) ; 
int linalg_solve (MultiMatrix *B, 
MultiMatrix *X, 
void *linalg_workspace , 
ErrorMsg error_message) ; 

The MultiMatrix structure is used to represent any matrix used by the evolvers, and it is 
defined in multimatrix.h. One notable exception from table 1 is LAPACK for large dense 
systems, which could easily be implemented using the above scheme. 
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